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Abstract 

Nonequilibrium dynamics of biomembranes with active inclusions is considered. The inclusions 
represent protein molecules which perform cyclic internal conformational motions driven by the 
energy brought with ATP ligands. As protein conformations cyclically change, this induces hydro- 
dynamical flows and also directly affects the local curvature of a membrane. On the other hand, 
variations in the local curvature of the membrane modify the transitions rates between confor- 
mational states in a protein, leading to a feedback in the considered system. Moreover, active 
inclusions can move diffusively through the membrane so that surface concentration varies. The 
kinetic description of this system is constructed and the stability of the uniform stationary state 
is analytically investigated. We show that, as the rate of supply of chemical energy is increased 
above a certain threshold, this uniform state becomes unstable and stationary or traveling waves 
spontaneously develop in the system. Such waves are accompanied by periodic spatial variation 
of membrane curvature and inclusion density. For typical parameter values, their characteristic 
wavelengths are of the order of hundreds of nanometers. For traveling waves, the characteristic 
frequency is of the order of a thousand Hz or less. 

PACS numbers: 87.16.Dg, 05.40.-a, 05.70.Np 
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I. INTRODUCTION 



Membranes play a fundamental role in various functions of living cells, providing spatial 
compartmentalization and being essential for signal transduction and inter- or intracellular 
transport While in vivo experiments with biomembranes are difficult because of their 
strong coupling to cytoskeleton and cytoplasm, investigations with synthetic membranes and 
vesicles formed by lipid bilayers can be performed [2|. Equilibrium morphology transitions 
in vesicles are well described by the elastic membrane theory [31]. However, it becomes in- 
creasingly clear that, under standard physiological conditions in a living cell, biomembranes 
are far from the state of thermal equilibrium. 

Lipid bilayers that constitute a biomembrane are usually including a large number of 
protein molecules which may represent ion pumps, enzymes or other molecular machines. 
The common property of active protein inclusions is that they undergo cyclic conformational 
changes, with each next cycle initiated by binding of an energy-bringing ligand, typically an 
ATP molecule. Active conformational motions in the inclusions are coupled to membrane 
dyua,„,cs. thus bringing a membrane away from the state of thermodynamie equUibrium^ 

Theoretical investigations of biomembranes including active ion pumps have shown [4,] 
that their fluctuational power spectra must be very different from that of the equilibrium 
membrane systems. It has been experimentally demonstrated that active inclusions renor- 
malize bending rigidity [s], ^ and surface tension [7[ of biological membranes. In several 
theoretical studies, it has been furthermore predicted j^, 9] that instabilities of membrane 
shapes may be induced by the nonequilibrium activity of the inclusions. Such instabilities 
could not however so far been observed in the experiments. 

The original theoretical analysis of active membrane instabilities has been performed 
in the framework of a general phenomenological approach where the details of operation of 
individual protein inclusions did not play a role. Effects of conformational transitions in 
active protein inclusions have been explicitly considered in the previous publication j^, [l^ 
where the equations for membrane dynamics have been complemented by a kinetic equation 
for conformational transitions in molecular inclusions. Only inclusions with two different 
conformational states have been previously studied. 

However, as we point it out now, models with two-state inclusions cannot describe generic 
properties of nonequilibrium membrane systems. Indeed, when only two discrete states are 
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present in active element and stochastic transitions between them take place, such transi- 
tions always satisfy the condition of detailed balance. Therefore, for such elements, there 
is no principal difference between the states of thermal equilibrium and nonequilibrium 
steady states, the latter being always describable as equilibrium states with some effective 
temperature. 

In order to see true differences in the behavior of equilibrium and nonequilibrium mem- 
branes, inclusions with at least three internal conformational states must therefore be con- 
sidered. For such inclusions, cyclic probabilities flows indicating the absence of the detailed 
balance are possible. The intensity of such flows provides direct characterization of the 
degree of deviation from thermal equilibrium. 

Having this in mind, we construct here a complete kinetic description for the membranes 
with multiple-state inclusions, representing a generalization of the previous theories j4,[5|, Q]. 




For the general case of the -ft"-state inclusions, a detailed stability analysis is analytically 
performed. In addition to the long-wavelength static instability due to the negative effective 
surface tension, which has previously been known, two new instabilities are found. They are 
the Turing-type static instability with a finite spatial wavelength and the oscillatory wave- 
type instability with a finite wavelength. The characteristic wavelengths of both instabilities 
are of the order of hundreds of nanometers; the characteristic frequency of the patterns 
emerging as a result of the wave instability is up to a thousand Hz. 

Below in Section II the theoretical model is formulated. First the expressions for the free 
energy of the system are given and then the dynamical evolution equations are constructed. 
In Section HI, the slow-time limit of the model, valid when the time scales of intramolecular 
processes are much shorter than the time scales of the merging membrane patterns, is 
introduced and discussed. The stability analysis of the uniform stationary state is performed 
in Section IV. The results are finally discussed in Section V. 

II. THE MODEL 

A. Free energy and the equilibrium state 

We consider a membrane composed of lipids and active inclusions with K internal states. 
In this study, we assume that all inclusions have the same orientation with respect to the 
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membrane. The situation where inclusions with both up and down orientations are present 
can be considered in a similar way. Before discussing the dynamics, we first consider the 
free energy of this system. At coarse-grained level, its free energy is given by the sum of 
the elastic energy of a lipid bilayer, the free energy associated with the inclusions, and the 
energy of coupling between the inclusions and the elastic membrane, 

F Fmem ~l~ Pine ~l~ Fcoup- (l) 

For a quasi-flat membrane, the elastic free energy of a lipid bilayer is 

Fmem =\j + K{V'hf + h{V'hf] , (2) 

where h is the height of the membrane relative to the reference plane, r = xi + is the 
two-dimensional position vector in the reference plane, and V = dji + dy'^. The coefficient 7 
characterizes surface tension of the membrane, k specifies the bending rigidity of the mem- 
brane, and the coefficient ^2 characterizes its rigidity with respect to the spatial variation 
of the curvature. Below, the last term is taken into account only for the short-wavelength 
instabilities, where it becomes essential and determines the instability wavelength. Since the 
system under consideration is a biomimetic membrane, we take the surface tension of typical 



cell membranes as the upper limit for the magnitude of 7, i.e, 7 ^ 5 x IQ^^ksT /nm^ 
The coefficient k typically is k ~ lOkeT [l^, and ~ c^^^ ~ 300 keT ■ nm^ (a ~ 5nm is 
the linear size of an inclusion). 

The free energy associated with the inclusions is 



no / d^r I ksT 



^ In (/)^ + (1 - ^ 0^) ln(l - ' 



where < 0q, < 1 is the dimensionless surface density of inclusions in the state a {a = 
1,2,3, ...,K), Ea is the energy of an inclusion in the state a, no = In this paper, we 

consider the case of the low inclusion density and neglect therefore possible direct lateral 
interactions between the inclusions. 

The coupling between the inclusion density and the local membrane curvature is 

Fcoup= j (f rKYca(j)aV^h, (4) 

where Cq are the coefficients specifying the strength of coupling between inclusions in different 
internal states a and the local membrane curvature. When Cq, > 0, membrane regions rich 



in the inclusions in the state a tend to have V^/i < 0. The coefficient Cq has dimension 
of inverse length, ~ ASa/EcZa where Sq, is the average cross-sectional area, AS^ is the 
difference between the outer and inner leaflet surface area of a state-a inclusion , and 
is the thickness (in the z direction) of the inclusion (see Appendix for the explanations). 
Typically, is of the order of < O.lnm^^. 

The equilibrium state of the system satisfies the conditions 5F/5h = and 6F/6(j)a = 0. 
This leads to the equations determining the equilibrium state, 

KV^/i-7V2/i-A;2V^/i + fi:^c«VVa = 0, (5) 

a 

and 

" +noK = 0. (6) 



Eq. ([6]) gives us 



In 0a - ln(l - y^0/3) 





0a = (1 - J] 0/3) exp [-{E^ + n^^KCaN''h)/kBT] . (7) 

/3 

Therefore the local total density of the inclusion in thermal equilibrium is simply 
, ^ , ^ Eg exp [-{E^ + nQ^KCc,V'^h)/kBT] 
V l + Y.c.^wHEc. + n^'t^c^V^h)lkBT\ 

Note that the total number of inclusions present in the membrane is N = no J (Pr (pf It 
depends on E^, implying that Ea includes the contribution from the free energy difference 
between an a-state inclusion in the membrane and an a-state inclusion in the solvent. 

Already the equilibrium state of the system may exhibit complex structural pha ses such 
as the undulated lamellar phase, the hexatic arrangement of inclusions, etc (see [iSj]). These 
phases are possible because of the strong coupling between the inclusion density and the 
local membrane curvature, and it is has been argued that they are important for biological 
membranes [l^. Because we are primarily interested in the nonequilibrium behavior, we 
limit however our study below to the situations where the uniform flat state of the membrane 
is stable under equilibrium conditions. 



B. Evolution equations 



If a membrane is impermeable, it can only move together with the solvent flow. A 
(partially) permeable membrane can be dragged through the solvent. Thus, the equation 
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for the temporal evolution of h{r, t) is 

' 5F 



dth = V ■ fi — A 



p 



a, /3 



(9) 



where v(r) is the local solvent flow and n(r) is the unit normal vector of the membrane 
surface. Because of the permeation, the solvent can "leak" from one side of the membrane 
to the other. If the force / is locally applied at a permeable membrane, it gets dragged 
through the solvent with the velocity v = Xpf, where Xp is the permeation constant. The 
local force / is given by the expression in the brackets in equation ([9]). In addition to the 
first contribution, directly coming from the free energy, it also includes the second term of 
purely kinetic origin. When an inclusion changes its internal conformational state, some 
force acting on the membrane is generated. Such average force corresponds to the second 
term inside the brackets. Here, fc^/? is the transition rate constant for an inclusion to go from 
state /3 to state a and P^^ is the momentum that an inclusion delivers to the membrane 
during a (3 —>■ a conformational change. Typically ranges from pN ■ fis to pN ■ ms. In our 
analysis, we shall neglect thermal noises acting on the membrane, which should be generally 
present in the membrane evolution equation. 

The local density of inclusions in the state a changes with time because of the lateral dif- 
fusion and the conformational transitions between the states, so that the evolution equations 
for the surface densities of inclusions in the state a are 

6F 

dt(t)a = ^aV^— + V [k^pcPp - kp^(f)^] , (10) 

where is the mobility coefficient for the state-a inclusions. 

Hydrodynamic flows in the solvent are described by the modified Stokes equation in the 
limit of low Reynolds number that takes into account the elastic stress due to the membrane 
and the forces resulting from the inclusion conformational changes. We have 

6F 

= -Vp -d,pz + r/(V^ + (9^)v - 5(z - h)—n 

oh 

+ Yl kap^oP.p<Pp [(1 + z^:l'^^h)5{z - z^:^p -h)-{l- ziyh)6{z + - h)] n, 

(11) 

where r] is the viscosity of the solvent. In this equation, 5F/5h corresponds to the pressure 
jump across the membrane due to the membrane elasticity. The last term comes from 
conformational changes of the inclusions. 



Because the dynamics of an inclusion in the fluid is strongly overdamped, an inclusion 
exerts a zero net force on the solvent during its conformational transition from the state 
j3 to the state a. The simplest description for the effect of an inclusion conformational 
change on the solvent flow is therefore provided by a model of a force dipole. z^^^ and zf^ 
in Eq. (ITT]) are characteristic lengths (both on the order of nm) for the force distribution 
during a conformational change of an inclusion from the state (3 to the state a. As has been 



previously pointed out in Ref. [l5|, but neglected in Refs. jol and 
of the inclusions lie in the surfaces parallel to the membrane (see Fig. [1]) the local density 
of the force centers in the plane h + z^^^ is (1 + z'^^pV'^h)5{z — z'^^ — h) and the local density 
of the force centers in the plane h — z^'^ is (1 — z^'^ V^/i)5(-2 + z^^^ — h). 
In addition, because the solvent is incompressible we have 

V ■Y + d,v, = 0. (12) 



Ml, because force centers 



C. Microscopic reversibility and detailed balance 

Microscopic reversibility imposes several constraints on kap, P^^, Pa/B, z'^p, and zf^. The 
local membrane curvature affects the energy landscape of the conformational states of an 
inclusion, thus in general the transition rates kap should depend on the local membrane 
curvature. To the order of V^/i we have 

kap = ^%k^!l{l + n%LpS/'h). (13) 

Here, the parameters k'^^'^ and lap correspond to passive inclusions, in the absence of external 
energy supply; k^^^ is the equilibrium transition rate and l^p represents the characteristic 
length associated with the curvature dependence of kf^. The parameters ^2^"^ and Vt^^^^ are 
introduced to describe changes in these properties when inclusions become active, i.e. when 
the energy is supplied to the inclusions. Thus, for passive inclusions Vt^^^^ = = 1 and 
the transition rates take the equilibrium form k^p = fc^^ = k^^pi^ + ^ap^'^h)- 

The detailed balance condition implies (to the order h) the following relationship between 
the equilibrium transition rates /c^^ and fc^^,, 

kpa e-(^/3+"o'''':/3V2/i)/fcBr \^ noksT J ' ^ ' 
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The second expression is valid as long as K,{ca — Cj3)'S7'^h/nokBT — [K,{ca — cp) /uokBT]/ R <^ 1, 
where R — | V^/i|~^ is the radius of the local membrane curvature. The above condition yields 



k 



,(0) 



U(0) 



(15) 



and 



a/3 ~ ''13a 



Iftn — 



noksT 

When the system is passive, there can be no active permeation. Therefore, we have 



(16) 



ad 



h.e pm 



0. 



(17) 



The requirement of the detailed balance k^gcpg — k%^(f)^^ yields 



(18) 



This relationship is a direct consequence of the microscopic reversibility of an a ^ transi- 
tion: the momentum transfer from an inclusion to the membrane during an a — /3 transition 
has the same magnitude but opposite direction as the momentum transfer from an inclusion 
to the membrane during a /3 ^ a transition. 

Two other conditions that must be satisfied for passive systems are 



in)' 

af3 



MY 

^al3 



)0/3 = 0. 



Again, because of the detailed balance, these two conditions lead to 



Pa/3 — —PjSa: 



(19) 



and 



Jn) _ (u) (d) _ (d) 



(20) 



These conditions are also direct consequences of microscopic reversibility. 
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D. Final evolution equations 



By using equations (fTT!) and (fT2l) and applying the Fourier transformation in the co- 
ordinate space, the flow velocity v can be expressed in terms of the hight variation h. 
Substituting the resulting expressions into ^ and (ITUi) . final evolution equations for the 
membrane hight and the inclusion density in different conformation states are obtained. 
They have the form 



6F 
Jh 



m 



a/3 



a/3 



^a/3 "T ^a/3 - 



a/3 



{uY 
^a/3 



^a/3 2 



(q,t), 

(211 



and 



5F 



dt(f)ai(l,t) = -Maq'^ — {ci,t) + ^ iK/3(l)l3 " A;/3„0c,) (q,t). 



(22) 



Here /i(q, t) and 0a (q, t) are the coefficients of the Fourier expansion of the local hight h{r, t) 
and the local densities 0a(r, t) over the plane waves exp(zq- r). 

The contributions to the membrane dynamics given by Eq. (1211) come (i) from the direct 
momentum transfer from the inclusions to the lipid bilayer, described by fca/3'^o-P(^0/3; (h) 
from the action of the active force dipoles due to the coupling between the local membrane 



"a/3 



+ ^i/])'?^^; and 



curvature and the force density distribution, described by kapn^Pc 
(iii) from the action of the active force quadrupoles due to the up-down asymmetry of the 
shape of the inclusions, described by kai3noPai3iz^ai3 ~ ^ap )Q'^0/3/2- 

Evolution equations ( l2Ti) and ( l22l) provide the final concise description of the membrane 
dynamics coupled to teh kinetic transitions in active multiple-state inclusions diffusing within 
the membrane. The physical meanings and typical magnitudes of the parameters entering 
into these equations and the respective expressions for the free energy are summarized in 
Table I. 

Equations (I^T]) and (12^ represent the generalization of the equations previously obtained 
in the theories sl, [l5| that neglected conformational changes inside the inclusions, where 
active momentum transfer from the inclusions to the membrane js] and active force dipoles 



and active force quadrupoles 



15| were both taken into account. 
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III. THE LIMIT OF SLOW MEMBRANE DYNAMICS 



Membrane motions with long wavelengths are slow. The characteristic time scales of such 
motions are much larger than the characteristic times corresponding to the the kinetics of 
internal transitions inside the inclusions. In this situation, which we consider further in our 
study, the steady nonequilibrium distribution over the internal states of active inclusions 
should adiabatically following the changes in the local membrane shape. 

For example, let us consider biomembranes with active ion pumps, such as bacteri- 
orhodopsin (BR) and Ca^'''-ATPase, which have used in the experiments fl, y, Q. The 
typical timescale for conformational transitions inside the cycle in these ion pumps is of 
the order of several milliseconds. The lateral diffusion constant for an an ion pump inside 
a membrane is on the order of 1/im^/s. Therefore, within a single active conformational 
cycle such an inclusion would move over the distance of about several tens of nanometers 
inside the membrane. Hence, if the membrane shape varies on a scale of several hundreds 
of nanometers or longer, the adiabatical approximation for the inclusion kinetics shall be 
applicable. In the steady state, local inclusion densities satisfy the equations 

"^{^(3(1)13 - ki3a(pa) = for all a. (23) 

In this paper, inclusions with cyclic transitions are considered (Fig. 2), so that k^p = for 
/3 — a 7^ ±1. Therefore, Eqs. (l23l) give us 

k2l4>l - ki2(t)2 = k324>2 - k23h = k^sh - ^34</>4 = = -^^s (/^a/3 )</>*• (24) 

Thus, conformational transitions are characterized by a steady conformational current Ks4>t- 
When Kg > the conformational current is along 1— >-2—>-3— >■...— >-K—>-l. The presence 
of such a flow (if Kg is not vanishing) indicates the absence of detailed balance and the 
deviation from the state of thermal equilibrium for the inclusions. 

Slow dynamics of the system is described by two variables: membrane height h and the 
total density of inclusions (pt = Xla ^a- From Eq. ( l22l) . we obtain 



9i0i(r, t) = J2 Ma^^ [kBTno In + Kc^V^h^ . 



(25) 



Here 0q, is a function of 0t and K^, we will discuss this function later after deriving the 
linearized equations. To find the equation for h, terms related to active forces are simplified 
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by introducing the active momentum transfer, the active force dipole and the active force 
quadrupole associated with an inclusion. From Eqs. f[T8l) . f|T9|) . and flMl) . we obtain 

ko.pP'^pC^p = {PTi + - + Pk k-1 + A"^)(fc2i0i - fci202) = Prir,(fc„^)0„ (26) 

0/3 



P2l(4? + 4f ) + ^32(4? + 4?) + - + PlKizfi + (^2101 - k 



12V2) 



PaKs 



(27) 



and 



MY 



p ^21 ^21 I p ^32 ^32 I I p ^li^ ^IX 
-'21 H -^32 V ■■■ + r\K 



(A;2101 - /Cl202) 



sVt- 



(28) 



Here, we have defined the active momentum transfer from an inclusion to the membrane 
during one cycle = P^ + P^ + ... + P[^, the net active force dipole of an inclusion during 



one cycle pa = P2i(4i^ + + -^32(42^ + 420 + ••• + Pixiz^'^x + ^ik)^ net active 



force quadrupole of an inclusion during one cycle Qa = P2i(4i'' ~ -2^21'' )/2 + P32(-2: 

^32 



MY 



MY 

32 



)/2 + ... + Pik{z^ix — z["i^ )/2. Note, that the respective terms are non-vanishing and 
thus contribute to the membrane dynamics only when the detailed balance is violated and 
the cyclic conformational current Ks(j)t is present. The equation of motion for /i(q, t) can be 
obtained by substituting Eq. ([26]), ([27D, and Eq. ([28]) into Eq. yielding 

'6F 



Jd) 



dth{ci, t) 



-X 



6h 
Jh 



+ noPTK 



sVt 



+ noPAKsCpthq"^ + noQAKsM"^ 



(q,t)- 



(29) 



Eq. (1291) indicates that the membrane acquires a nonzero velocity due to active perme- 
ation, as one can see from dth = —\pnQP^Ks(t)t at g = 0. For a permeable membrane this 
drift provides clear manifestation of non-equilibrium processes in the system under consid- 
eration. 
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IV. NONEQUILIBRIUM INSTABILITIES 



To analyze possible instabilities of the flat membrane, we should linearize the equations 
of motion by expanding h{r,t) and (f)t{'r,t) around their uniform solutions /io(^) and (pto, 

h{r,t) = ho{t) + 6h{r,t), 

0t(r,t) =0to + ¥t(r,t). (30) 

Note that Sh(c[,t) = h{q,t) and (5<^t(q, t) = </)t(q, t) for g 7^ 0. 

First, we consider the linear elasticity of the system. One can express 6F/6h as 

SF 



Sh{r) 

At small V^/i, we can further express (f)a/4>t as 

^ = ^^{l + l^^V'h + ...), (31) 

where tpa is dimensionless and l^a is a length that characterizes the curvature dependence 
of (pa/(pt- Because all coordinate dependence of 4>a/<Pt comes from the local curvature 
dependence of the transition rates, both ipa and l^a are independent of r, but they de- 
pend on Q^^p and f^^^^, i.e., on the strength of the external energy drive. From Eq. (!3T!) . 
V^(0t0„/0t) = ^/'aV^t + ipJ^a4>to'^'^h + .... As a resuh 



6F 

^ ' a a 

= KeffV^h--fV^h-k2V^h + Keffc''^ffV^(j)t. (32) 

Here we have introduced the effective membrane bending rigidity Hejf = k(1 + 
J2a (^a4^oJ"ii)a<t'to) and the effective inclusion-membrane elastic coupling constant Cgj^ = 
^a^oii'ai^/ i^eff- Note that, besides c^. Kg// and c^^j also depend on and fi^^^j through 
ifja and /^a. These effective elastic moduli are different from their equilibrium values because 
the relative inclusion population (pa/'Pt is different from the equilibrium distribution. 
Next we consider the equation for 5(j)t. From Eq. (I3T!) . 

2 A„ (t>a 



In — ^ = VM In — + In 



1 V 1- 
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Substituting the above expression into Eq. (1251) . one finds 

It is convenient to introduce the effective mobihty Mg// = J2a^a^ effective suscep- 
tibihty of the inclusions Xeff = ""-o^s^ [1/0*0 + 1/(1 — 0io)] and the effective inclusion- 
membrane elastic coupling constant c^jj associated with inclusion diffusion, MeffKeffC^^j = 

Mq, {nca + riokBTl^a)- Similar to n^ff and Cgjj, c^jj- depends on and ^^l- 
equation for 5f0i(q, t) can be written as 



dtMcL, t) = -M^ff (^XeffQ^M^, t) - '^effctffq'^h{ci, t)^ . (33) 

To obtain the linearized equation for dh{q,t)/dt, we expand the conformational current 
to the order of h as Kg^ka/s) = Kq{1 + l^'^'^h + ...) and obtain 

Ks{kap)(l)t = Ko(l)to + Ko6(f)t + KolK^to'^^h, (34) 

where Kq and Ik are both independent of r, but they depend on Q^^^ and fi^^j. Substituting 
Eq. ( |32l) and into Eq. (l29l) . the linearized equation of motion for /i(q, t) in the momentum 
space is found, 

dMq,t) = -\ [(7^g2 + ^ ^^^6) ^^q^^) _ ^^.^^^c^^q^ - rioPTKo) 0*(q,t)} 

~4^ [(7ag^ + /^ag"^ + ^^2?^) h{q, t) - KejfCaq^(j)t{q, t)] . (35) 

In this equation, following notations have been introduced: 

7^ = 7 - noPTKolK<Pm (36) 
is the surface tension renormalized by active forces that inclusions act on the membrane, 

7a = 7 + noPAKo(j)to (37) 
is the surface tension renormalized by active force dipoles, 

f^a = l^eff — noQAKolK(f)to (38) 
is the bending rigidity renormalized by curvature- induced active quadrupoles, and 

l^effCa = l^effCeff - UoQaKq (39) 
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defines the renormalized inclusion- curvature coupling due to curvature-independent active 
force quadrupoles. 

Since experiments are usually carried out for the membranes with negligible permeabil- 
ity , we consider below only instabilities of impermeable membranes with Ap = 0. 
It is convenient to write the derived the equations in the matrix form 



where 



Arj 

Mh^iq) = -^i^effCa, 

M^^{q) = -M.ffXeffq'. (41) 

The definitions and physical meanings of the parameters in Eq. fHOl) fj4T|) are summarized 
in Table II. 

To analyze the linear stability of the system, we seek solutions of Eq. (1401) of the form 
h{q,t), 0j(q, t) ~ e'^*. The characteristic equation for A is 

A' - {Mhh + M^^)X + {MhhM^^ - M^hMh^) = 0. (42) 
Let the two solutions of the above equation be Ai(q), A2(q), then 

Ai + A2 = Mhh + M^^ = -^ {-iaq + 4:r]MeffXeffq^ + Kaq^ + k2q^) , 

A1A2 = MhhM^^ - Mh^M^h = ^""^""^ [laq^ + /^ag' + k2q') , (43) 

where 

2 4> 

Ka = —. (44 

Xeff 

The following instabilities can occur in this system: 

(i) A long-wavelength instability at g = when one of the A becomes zero at g = and 
both Ai, A2 are negative for any nonzero q. 

(ii) An oscillatory wave instability (a Hopf bifurcation with a finite wavenumber) 
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at g = Qh, when Ai + A2 = Mhhiq) + M^^{q) < for all q except q = qh and 
A1A2 = Mhh{q)Mtf)tf){q) — Mh^{q)M^h{,<i) > for all q. The frequency of the wave at the 
onset of the instability is uj/2-k = a/ Ai(g/i)A2(q'/i)/27r. 

(iii) A static Turing-type instability when Ai + A2 < for all q and A1A2 > for all q except 
q = Qs- 

The typical magnitudes of the parameters that are important for stability analysis are 
summarized in Table III. Figure [3] shows the phase diagrams for systems with typical pa- 
rameters. Varying the parameters does not change the qualitative features of the phase dia- 
grams. The ordinates in Fig. [3] represent the strength of curvature- dependent pumping, and 
the abscissas represent the strength of curvature-independent pumping. The state of thermal 
equilibrium corresponds to the origin of coordinates; as the system is driven out of equilib- 
rium, the state of the system moves along a straight line with the slope Ik /a- Generally, 
coefficients ^e//, c^ff, and c'^jj- also vary as the system is being moved away of equilibrium; 
for simplicity we take however constant values Kg// = ^ksT, and c^jj = c'^jj = O.lnm^^ 
when the phase boundaries below in Fig. [3] are determined. 

A. The long-wavelength static instability 

According to Eq. (Il3il . this instability occurs when 7a = 7 + '^oPA-f^o^to < 0, i.e. when 
the effective surface tension becomes negative. As shown in Fig. [3l for systems with < 
this instability is possible at sufficiently large Kq. Using as typical values \pa\ ~ nm ■ pN-/is 
to nm ■ pN ■ ms, Kq ^ ms~^ and Uq = a^^ ~ l/25nm^, we find that to get 7a < for 
relatively low inclusion density, such as ipto ^ 10~^, a vesicle with the bare surface tension 
7 ^ l{)~^kBT /mv? is needed. In the experiments the measured surface tension of a 
vesicle with passive BRs was 7 ~ 10~^A;BT/nm^. Thus, it should be possible to observe this 
long-wavelength instability in the experiments. 

B. The finite-wavelength instabilities 

Besides the long-wavelength instability corresponding to negative 7a, there are also in- 
stabilities when Ha or ka become negative. If 7a is positive, these instabilities correspond 
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to the growth of modes with finite wavelengths, because long-wavelength fiuctuations are 
suppressed by the membrane tension. 

(i) When 

f^a = -2^/^^, (45) 
the static Turing-type instability takes place. Its characteristic wavenumber is 

Qs = ila/k^f'. (46) 
As follows from Eqs. (EZD, (JSHD and (|44]), Eq. (|45]) is equivalent to 

/ ^ \ ^2 ^(l> 

UoQaKo lK<Pto - - 2^(7 + rioPAKMh = /^e// - ■ (47) 

V Xeff J Xeff 

As seen in Fig. HI two conditions are essential for the occurence of this instability. The 
first is that the difference Keff — K^jj-cfjc'^jj-fxeff should be small, implying that the curved 
membrane regions are strongly attracting the inclusions. The second property is that the 
combination noQAKQ{lK(pto — f^effC^ff/Xeff) should be large and positive, implying strong 
curvature-induced active force quadrupoles. In the phase diagram in Fig. [3], the Turing-type 
instability is located where KqIk is large. 

The characteristic wavelength at the onset of the Turing-type instability can be estimated. 
Taking 7a ~ 7 ^ 5 x lO^'^/^^T/nm^, ^2 ~ SOOksT ■ nm^, we find that this finite-wavelength 
instability has a characteristic wavelength of 27c /qs = 271/(7^/^2)"'^''^ ~ lOnm — lO^nm. 
Numerically computed wavelengths for this instability are shown in Fig. [51 Note also that, 
according to Fig. [3] and Fig. [5l for < only the membranes with inclusions that have 
large Iq (i.e., the activity of these inclusions is induced by membrane curvature) may exhibit 
the Turing-type instability, otherwise long- wavelength instability is taking place. 

(ii) When = k^, the oscillatory wave instability with the characteristic wavevector qh 
takes place. The values of and qh are given by the solutions of two equations. 



Qh 

and 



r \ 2 



24 + + ^_^MfMl = 0, (49) 
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Because qh > 0, Eq. (HHi) requires to be negative. Generally, the above equations for 
and Qh have to be solved numerically. Fig. [3] shows that this instability may occur only 
when Pa, Qa, and c^^j, cf^j all have the same sign, and the curvature-induced pumping is 
not too strong (i.e., Ik is not very large). This condition suggests the following picture for 
the mechanism of wave instability. As shown in Fig. [6], first, inclusions are attracted to the 
regions with their preferred membrane curvature; but the active forces from the inclusions 
in these regions produce curvature that the inclusion dislike (thus Ik cannot be large). 
As a result, the inclusions move back and forth, chasing regions with preferred membrane 
curvature. This leads to an oscillatory instability with the characteristic wavelength which 
is large as compared to the molecular scale. 

The characteristic wavelength at the onset of wave instability is about 271 /qh ~ lOOnm. 



Substituting to Uh = a/ Ai(gh)A2(q'h), the characteristic time at the onset of wave insta- 
bility is estimated to be of the order of a millisecond (see Fig. [7] for numerically computed 
characteristic time at the onset of wave instability). This is already on the border of the 
validity of the approximation of slow membrane dynamics, because such characteristic time 
is comparable with the time scale of conformational transitions in individual active inclu- 
sions. Therefore, our analysis of the wave instability should be viewed as only providing 
an indication that such an instability may exist. The accurate quantitative analysis of this 
instability should be performed without assuming the slowness of membrane dynamics; it 



will be undertaken in a separate publication. [l6|] 



V. DISCUSSION 



We have constructed a complete kinetic description for membranes with multiple-state 



active inclusions. This work generalizes previous theories on active membranes [8|, 19|, Il5 |. 
Although we assume that all inclusions have the same orientation with respect to the mem- 
brane, the situation where inclusions with both orientations are present is a simple extension 
of this work. The active forcing resulting from the conformational changes of the inclusions 
affecting membrane dynamics through a direct force acting on the membrane (active per- 
meation), a force dipole due to density asymmetry induced by membrane curvature, and a 
force quadrupole, are all present in Eq. ( l2Ti) . Comparing to previous theoretical works on 



active membranes, Eq. ( ]2T]) directly connects the active forcing with each conformational 
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change of the inclusions. 

In the hmit of slow membrane dynamics, the distribution over the internal states of ac- 
tive inclusions adiabatically following the changes in the local membrane shape. Our theory 
recovers the active force dipoles and active force quadrupoles introduced in previous theo- 
ries js], [isl. Furthermore, the linearized equations of motion (HQ]) (141 p reveal that the effective 
elastic moduli of the membrane are renormalized both by the active forcing resulting from 
inclusion conformation changes and the distribution over the internal states of the inclusions. 
Thus, our linearized equations of motion can be applied to analyze the experimentally ob- 
served renormalization of the membrane rigidity jo] and surface tension [3] by inclusion 
activities. It would be interesting to see how the pumping effect due to the nonequilibrium 
forcing and the non-pumping effect due to the nonequilibrium inclusion distribution over the 
internal states contribute to the observed mechanical properties of an active membrane. |l^ 

We also show that when the system is driven sufficiently far away from the equilibrium 
state, a uniform membrane becomes unstable and, besides a long-wavelength instability 
arising from negative effective surface tension, static Turing-type or traveling waves spon- 
taneously develop in the system. Figure [3] shows that the possible instabilities for an active 
membrane depend on the signs of the active force dipole pa, active force quadrupole Qa, and 
whether inclusion activity is sensitive to the local membrane curvature (described by the 
magnitude of the parameter Ik)- For inclusions whose activity is highly sensitive to the local 
membrane curvature, it is possible for a static Turing-type instability with a characteristic 
length of the order of hundreds nanometers to take place. For inclusions whose pa, Qa, 
renormalized elastic constants c^jj and c^j^j all having the same sign and activity is not very 
sensitive to the local membrane curvature, it is possible for a wave-instability to take place. 
The characteristic length of the traveling wave is of the order of hundreds of nanometers 
and characteristic frequency is of the order of a thousand Hz or less. 

For traveling waves with frequency of the order of a thousand Hz, the characteristic 
time is already comparable with the time scale of conformational transitions in individual 
active inclusions. Therefore, in the accurate quantitative analysis of this type of waves, the 
inclusion distribution over the internal states should not adiabatically follow the changes 
in the local membrane shape. Such analysis should start from Eq. (12T]) . waves with high 
frequency and wavelength of the order of hundreds of nanometers should manifest themselves 
as the synchronization of inclusion conformational cycles, [l^ 
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We believe our theory will inspire new experiments on the active membranes. Once pa 
and Qa of the inclusions are identified through the measurements of the surface tension and 
the bending modulus of the passive and active membranes, the criterions for the instabilities 
can be deduced from our theory. The static Turing instability may be observed by optical 
microscope, and both stationary and traveling waves can be detected by dynamics scattering 
experiments. In the future, it is also possible to design "active vesicles" whose morphology 
and mobility are controlled by the activities of the inclusions. Such active vesicles may 
harvest the energy supply in the environment, turn chemical energy into mechanical energy, 
and act as membrane machines. 
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Appendix 

In this appendix we give a brief calculation that relates the shape of a conical state-a 
inclusion to its coupling constant Cq,. 

As shown in figure [H the outer radius, inner radius, and the thickness of a state-a 
inclusion are r°"*, r^", and respectively. The preferred radius of curvature Ra is related 
to these quantities by 

out in 
a a 



Ra + ^a/2 Ra — L/2 



That is. 



^ L(rr + C)/2 (rr + C)^^a/2 ^gla .gQ. 

where Sq, is the average cross-sectional area of a state-a inclusion, ASq, is the difference 
of cross-sectional area between the outer and inner surfaces of a state-a inclusion. The 
coupling constant Ca is simply Cq ~ l/Ra- From Eq. (I2])(jlj), regions rich in > inclusions 
tend to have V^/i < 0, this is illustrated in Fig. [91 
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TABLE I: notations in Eg. (|2D (|3D (ID (|2TD (|22]) 



symbol 


physical meaning 


typical value 


K 


membrane bending modulus 


~ WksT 


7 


membrane surface tension 


< 5 X IQ-^kBT/nm'^ 


Cq, 


inclusion-curvature coupling constant 


^ O.lnm^^ 


k2 


higher bending modulus for the membrane 


~ SOOfcer • nm^ 


a 


lateral linear size of an inclusion 


5nm 


no 


l/a2 


~ 0.04nm"^ 


V 


solvent viscosity 


~ 0.01 - O.lg/cm • s 


kaP 


inclusion conformational transition rate 


~ ms~^ — fis~^ 




momentum transfer to the membrane for /3 — > q inclusion 
conformational change 


~ pN • fis — pN • ms 




momentum transfer to the solvent for /3 — > a inclusion 
conformational change 


~ pN • yus — pN • ms 




characteristic lengths associated with 
active force during (3 ^ a transition 


~ nm 
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TABLE II: notations in Eq. (gOD 



symbols 


definition 


physical meaning 






represent the relative density of state-a inclusions 






effective bending rigidity in 6F/5h 






inclusion-curvature coupling in 6F/5h 




XIq, MainCa + nokBTl^a)/MeffKeff 


inclusion-curvature coupling in d(f>t/dt 






KQ(f)t is the curvature- independent conformational 
current 


Ik 




Kolkiy'^h)(j)to is the curvature-dependent 
conformational current 


la 


7 + noPAKo(l)to 


surface tension renormalized by active force dipoles 




l^eff - noQAKolK4>tO 


bending rigidity renormalized by curvature-induced 
active quadrupoles 


Ca 


{Keffdlff - noQAKo)/Keff 


inclusions-curvature coupling renormalized by 
curvature-independent active force quadrupoles 
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TABLE III: typical magnitudes of parameters that are important for stability analysis 



symbol 


physical meaning 


typical magnitude 


VA 


active force dipole 


\pA\/a < pN • ms 


Qa 


active force quadrupole 


IQaI/o^ ^pN-ms 


Ko 


characteristic time scale for curvature-independent 
inclusion conformational-change cycle 


|i^o| ~ ms 




characteristic time scale of curvature-induced 
inclusion conformational-change cycle 


\KolKV'^h\ < ms 


<t>to 


average areal fraction of the inclusions 


0to ^ 0.1 


Xeff 


effective susceptibility of the inclusions 


Xeff ~ kBT/mr? 




effective membrane curvature elastic modulus 






effective inclusion-curvature elastic coupling constants 





24 



Figure Captions 

• Figure 1. Force centers of the active force dipoles for /3 — >■ o; transitions are distributed 
on surfaces described hy h + z^^^^ and h — zf^^. These surfaces and h{x) are parallel 
surfaces, the areas of parallel surface elements from top to bottom are (1 — 2;^^^^ V^/i)(ZA, 
dA, and (l + zjfjv^/i)^^. 

• Figure 2. Schematics of the conformational transitions of an inclusion with five inter- 
nal states. When one considers the dynamics of the system at lengths large compared 
to 0(10 nm), the strength of the conformational current of an inclusion with charac- 
teristic time scale comparable to typical ion pumps is described Kg^t- 

• Figure 3. Phase diagrams for membranes containing active inclusions. 7 = 
10-^A;Br/nm2, hi = 250A;sr • nm^, = 0.1, M^ffXeff = l/^m^/s, 77 = lO^^^kg/m • s. 
Long-dashed curve indicates the onset of wave instability, short-dashed curve in- 
dicate the onset of Turing instability, solid line indicates the onset of long wave- 
length instability. The phase boundaries are determined by taking /tg// = S/c^T 
and Cg^j — c^^j: — O.lnm"^ at the onset of all instabilities, (a) Pa/ci — IpN • ms, 
Qa/cl^ — IpN • ms. For inclusions with Ir/cl greater than the slope of the dashed 
line, the system has Turing-type instability at large energy input, otherwise the sys- 
tem has wave instability, (b) Pa/o — IpN • ms, Qa/o>^ — — IpN • ms. In this case, 
only Turing-type instability occurs at high energy input, (c) pA/a = — IpN • ms, 
QaIo^ = IpN ■ ms. For inclusions with Ik/cl greater than the slope of the dashed line, 
the system has Turing- type instability at large energy input, otherwise the system 
has long wavelength instability, (d) pa/cl — — IpN • ms, Qa/cl^ — — IpN • ms. For 
inclusions with \Ik/(i\ greater than the absolute value of the slope of the dashed line, 
the system has Turing-type instability at large energy input, otherwise the system has 
long wavelength instability. 

• Figure 4. Mechanism for Turing-type instabihty: (i) inclusions are attracted to regions 
with preferred curvature, (ii) curvature-dependent active force quadrupoles generate 
membrane curvature that inclusions prefer. This positive feedback leads to a Turing- 
type instability. 

• Figure 5. characteristic wavelengths at the onset of instabihties for 7 = lO^^/csT/nm^, 
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k2 = 250/cbT • nm^ 0to = 0.1, MeffXeff = lAtm^/s, 77 = 10~^kg/m • s. Kg// = 
BknT, and c^^j = c^j^ = O.lnm"^. The characteristic wavelengths are on the order of 
lOOnm for all cases except at the vicinity of long-wavelength instabihties, where the 
characteristic lengths increase abruptly, (a) thin curves: p^/a = IpN ■ ms, Qa/o-^ = 
IpN • ms (dashed curve: wave instability, solid curve: Turing-type instability); thick 
curves: pa/o — O.lpN • ms, Qa/o? — IpN • ms (dashed curve: wave instability, solid 
curve: Turing-type instability), (b) dashed curve: Pa/o. — — IpN • ms, Qa/o? — 
IpN • ms; and solid curve: Pa/cl — — IpN • ms, Qa/cl^ — — IpN • ms. 

Figure 6. Mechanism for wave instability: (i) inclusions are attracted to regions with 
their preferred curvature (ii) in inclusion-rich domains, active force quadrupoles pro- 
duce membrane curvature that inclusions dislike, the inclusions leave inclusion-rich 
domains and the membrane configuration goes back to (i). 

Figure 7. Characteristic time scale at the onset of wave instabihty for 7 = 
IQ-^kBT/n-m^, = 250A;sr • nm^, = 0.1, M^ffXeff = l/^m^/s, 77 = IQ-^kg/m • s. 
Keff = 5/cbT, and c^jj = c^^j = O.lnm"^. Dashed curve: Pa/ci = IpN • ms, 
QaI^^ = IpN ■ ms; solid curve: pA/a = O.lpN ■ ms, Qa/cP' = IpN ■ ms. The magni- 
tude oHk/o, for the solid curve is indicated at the top of the frame, and the magnitude 
of Ik /a for the dashed curve is indicated at the bottom. 

Figure 8. A conical inclusion in state a, la is the thickness of the inclusion, r°"* is the 
outer radius, r^" is the inner radius of the inclusion. 

Figure 9. Regions rich in positive inclusions tend to have V^/i < 0. 
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